****************************************************************************
* File-Nale: 		itemanalysis.do
* Date:		 04/06/2020
* Author: 		Fred Batista
* Purpose: 		Analysis of bias in item difficulty
* Data used: 		lapop10irt2.dta
* Data Output:	None	*/
****************************************************************************

* setting edvreg

program define edvreg
version 6.0
tempvar svar dvv varCorr
syntax varlist [if] [in], Dvste(string) [Prop] [borjas] [cluster(varlist)] [hc3]   [Level(integer 95)]
local indvars "`varlist'"
if "`cluster'"~="" {
  local cluster="cluster(`cluster')"
}
gettoken depvar indvars : indvars
qui gen `dvv' = `dvste'*`dvste' `if'`in'
disp
disp in green "Sampled dependent variable regression"
disp
if ("`prop'" == "") {
      qui gen `svar' = sum(`dvv')
      local swsq = `svar'[_N]
      qui reg `depvar' `indvars' `if'`in', level(`level') `cluster' `hc3' 
      local dof = e(df_r)
      local sumsr = `dof'*e(rmse)*e(rmse)
      qui matrix accum XpX  = `indvars' `if' `in'
      qui matrix accum XpGX = `indvars' [iw=`dvv'] `if' `in'
      matrix tr = trace(inv(XpX)*XpGX)

      *following line modified by E. leoni
      local sigma2 = (`sumsr'-`swsq') / `dof'      
      if ("`borjas'"=="") {
        local sigma2 = (`sumsr'-`swsq'+tr[1,1]) / `dof' 
      }
      if (`sigma2' < 0) {local sigma2 = 0}
      disp in green "Estimated Sigma     = " in yellow %7.5f sqrt(`sigma2')
      disp in green "Mean Omega2         = " in yellow %7.5f sqrt(`swsq'/`dof')
      qui gen `varCorr' = (1/(`sigma2'+`dvv')) `if'`in'
        regress `depvar' `indvars'  `if'`in' [aw=`varCorr'],  `cluster' `hc3' level(`level')
    }
   else {
     tempvar e v
      qui reg `depvar' `indvars'  `if'`in', level(`level') `cluster' `hc3'
      predict `e' if e(sample), resid
      replace `e' = `e'*`e'
      qui reg `e' `dvv' `if'`in',  level(`level') `cluster' `hc3'
      if (_b[_cons]>0) {
         disp in green "Coeff. on variance from sampling = " _continue
         disp in yellow %7.4f _b[`dvv'] " (" %7.4f _se[`dvv'] ")"
         disp in green "   Estimated variance of epsilon = " _continue
         disp in yellow %7.4f _b[_cons] " (" %7.4f _se[_cons] ")"
      }
      else {
        qui reg `e' `dvv' `if'`in',   nocons level(`level') `cluster' `hc3'
         disp in green "Coeff. on variance from sampling = " _continue
         disp in yellow %7.4f _b[`dvv'] " (" %7.4f _se[`dvv'] ")"
         disp in green "             Variance of epsilon = " _continue
         disp in yellow %7.4f 0 in blue " (constraint imposed)"
       }
      predict `v' if e(sample)
     regress `depvar' `indvars'  `if'`in' [aw=1/`v'],  `cluster' `hc3' level(`level')
   }
end


* renaming variables

* findit dm88_1

rename V1 provdiff

rename V2 termdiff 

rename V3 uspresdiff

rename V4 provdifferror

rename V5 termdifferror

rename V6 uspresdifferror

rename V7 numbprovince

rename V8 roundnumber

rename V9 lastchangeprov

rename V10 federal

rename V11 reelection

rename V12 fivesixyears

rename V13 parliamentary

rename V14 lastelection

rename V15 tourism09

rename V16 mercha09

rename V17 migstock10

** analysies

* nuimber of provinces

reg provdiff numbprovince lastchangeprov roundnumber

reg provdiff numbprovince lastchangeprov roundnumber, robust

reg provdiff numbprovince lastchangeprov roundnumber, vce(hc3)

reg provdiff numbprovince lastchangeprov roundnumber [weight=provdifferror]

edvreg provdiff numbprovince lastchangeprov roundnumber, dvste(provdifferror)

* length of term

reg termdiff reelection fivesixyears lastelection

reg termdiff reelection fivesixyears lastelection, robust

reg termdiff reelection fivesixyears lastelection, vce(hc3)

reg termdiff reelection fivesixyears lastelection [weight=termdifferror]

edvreg termdiff reelection fivesixyears lastelection, dvste(termdifferror)

* president of us

reg uspresdiff tourism09 mercha09

reg uspresdiff tourism09 mercha09, robust

reg uspresdiff tourism09 mercha09, vce(hc3)

reg uspresdiff tourism09 mercha09 [weight=uspresdifferror]

edvreg uspresdiff tourism09 mercha09, dvste(uspresdifferror)



**** Dif tests (critical t = 1.96)

gen id = _n

* provdiff

edvreg provdiff, dvste(provdifferror)

gen provt = (provdiff - 0.0381305)/(provdifferror + .2673002)

list id provt

gen provnorm = normal(provt)

replace provnorm = 1 - provnorm if provnorm > 0.5

list provnorm

rep

* termdiff 

edvreg termdiff, dvste(termdifferror)

gen termt = (termdiff -1*(-1.137161))/(termdifferror + .1306862)

list id termt

gen termnorm = normal(termt)

replace termnorm = 1 - termnorm if termnorm > 0.5

list termnorm

* uspres

edvreg uspresdiff , dvste(uspresdifferror)

gen usprest = (uspresdiff -1*(-.8383463))/(uspresdifferror + .114543)

list id usprest

gen uspresnorm = normal(usprest)

replace uspresnorm = 1 - uspresnorm if uspresnorm > 0.5

list uspresnorm

* see excel spreadsheet

